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Abstract.  This  article  presents  the  theoretical  motivation,  implementation  approach,  and  example  validation  results 
for  a  computationally  efficient  plume  simulation  model,  designed  to  replicate  both  the  short-term  time  signature  and 
long-term  exposure  statistics  of  a  chemical  plume  evolving  in  a  turbulent  flow.  Within  the  resulting  plume,  the  odor 
concentration  is  intermittent  with  rapidly  changing  spatial  gradient.  The  model  includes  a  wind  field  defined  over 
the  region  of  interest  that  is  continuous,  but  which  varies  with  location  and  time  in  both  magnitude  and  direction. 
The  plume  shape  takes  a  time  varying  sinuous  form  that  is  determined  by  the  integrated  effect  of  the  wind  field. 
Simulated  and  field  data  are  compared.  The  motivation  for  the  development  of  such  a  simulation  model  was  the  desire 
to  evaluate  various  strategies  for  tracing  odor  plumes  to  their  source,  under  identical  conditions.  The  performance  of 
such  strategies  depends  in  part  on  the  instantaneous  response  of  target  receptors;  therefore,  the  sequence  of  events 
is  of  considerable  consequence  and  individual  exemplar  plume  realizations  are  required.  Due  to  the  high  number  of 
required  simulations,  computational  efficiency  was  critically  important. 
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1.  Introduction 


Olfactory-bcLsed  mechanisms  have  been  used  to  explain  a  variety  of  biological  orientation  behaviors 
for  resource  location  [1,  2,  3,  4]:  homing  by  Pacific  salmon  [5]  and  green  sea  turtles  [6];  foraging 
by  Antarctic  procellariiform  seabirds  [7],  lobsters  [8,  9,  10],  and  blue  crabs  [11];  and  mate  finding 
and  foraging  by  insects  [12,  13,  14].  Typically,  olfactory-based  mechanisms  proposed  for  biological 
entities  combine  orientation  behaviors  at  long-range  based  in  part  on  olfaction  with  a  local  search 
in  the  vicinity  of  the  source  governed  by  multiple  sensor  modalities.  The  olfactory  portion  of  the 
search  strategy  is  contingent  on  the  fact  that  a  turbulent  flow  can  transport  above  threshold  parcels 
of  odor  over  relatively  long  distances,  making  feasible  the  detection  of  an  upflow  odor  emitter  over 
considerable  distances.  A  further  advantage  to  the  searcher  is  that  the  odor  signal  can  be  highly 
‘specific,  resulting  in  rare  false  positive  responses  [15].  The  flight  of  a  male  moth  along  a  wind- 
dispersed  plume  of  pheromone,  emitted  either  by  a  female  or  a  dispenser  charged  with  synthetic 
pheromone  is  particularly  well-documented  [12,  16,  17,  18,  19,  20].  Models  of  moth  navigation 
generally  propose  that  males  fly  upwind  when  they  sense  an  above- threshold  concentration  of 
pheromone.  This  mechanism  is  termed  ‘positive  anemotaxis.’  Because  in  a  rapidly-flying  insect, 
navigational  decisions  require  continual  review  of  the  sensor  by  the  insect,  the  fine-scale  plume 
structure  is  of  considerable  consequence. 

The  fact  that  various  biological  entities  use  olfactory  based  search  with  high  degrees  of  success 
has  prompted  research^  into  the  design  of  autonomous  vehicles  capable  of  olfactory  based  search 
and  chemical  plume  tracing.  Such  autonomous  vehicle  capabilities  have  applicability  in  searching  for 
environmentally  interesting  phenomena,  hazardous  chemicals,  and  pollutants.  In  addition  to  vehicle 
development,  sensor  development,  and  analysis  of  the  fluid  environment,  design  and  optimization  of 
the  chemical  plume  tracing  (CPT)  strategies  are  of  critical  importance  to  the  solution  of  this  prob¬ 
lem.  Evaluation  of  CPT  strategies  requires  hardware  and  software  evaluation  platforms.  During  the 
initial  design  stages,  software  evaluation  is  preferred,  since  such  tools  allow  competing  strategies  to 
be  evaluated  under  identical  conditions  for  various  environmental  scenarios.  The  tradeoffs  between 
field  testing,  wind  tunnel  testing,  and  simulation  testing  are  discussed  in,  for  example,  [33,  34]. 
To  calculate  reliable  performance  statistics,  numerous^  batches  of  such  evaluations  are  required. 
Therefore,  the  computational  efficiency  of  the  software  simulation  tools  is  a  key  concern.  To  ensure 
that  strategies  that  perform  well  in  simulation  will  also  perform  well  in  hardware  experiments,  the 


simulation  must  contain  the  key  features  that  complicate  the  problem  of  CPT.  This  article  presents 
a  plume  simulation  implementation  that  addresses  the  tradeoffs  between  computational  efficiency 
and  inclusion  of  realistic  complicating  features. 

Plume  models  can  represent  concentration  by  its  mean  or  probability  density  function.  For 
example,  dispersion  in  a  turbulent  medium  is  frequently  represented  as  a  Gaussian  distribution 
such  as  [21,  22,  23,  24] 

C{x,y,z)  =  +  ^)) 

where  C{x^y^z)  represents  the  average  concentration  profile  as  a  function  of  position.  In  this 
expression,  and  throughout  this  article,  x,  y,  and  z  denote  the  downwind,  crosswind,  and  vertical 
position  coordinates  relative  to  the  odor  source  with  x  positive  along  the  mean  wind  direction. 
Sy  and  Sz  are  the  standard  deviations  of  the  time-averaged  plume  concentration  in  the  crosswind 
and  vertical  directions.  Both  parameters  are  functions  of  the  downwind  position.  The  parameters 
of  the  functional  fit  may  be  based  on  theoretical  considerations  (i.e.,  Sy  =  and 

Sz  —  [23,  24]  or  selected  to  match  experimental  data  for  the  atmospheric  conditions 

of  interest  [22,  25,  26].  The  application  of  Gaussian  plume  concentration  models  is  well  established 
and  the  results  they  produce  are  well  understood.  Predictions  based  on  Gaiissian  plume  studies  yield 
results  that  match  experimental  results  reasonably  well  when  long-term  exposure  is  of  interest  and 
the  effect  of  the  exposure  is  a  linear  function  of  the  time  averaged  concentration  [27].  Probability 
density  plume  models,  however,  do  not  set  out  to  provide  information  related  to  the  short  time-scale 
signature  of  the  concentration,  the  instantaneous  peak  concentration,  or  the  intermittent  nature  of 
the  plume.  Although  the  fine-scale  characteristics  of  a  plume  are  of  little  interest  in  many  cases,  for 
example  in  some  kinds  of  regulatory  work  related  to  permitting  or  safety,  they  are  of  considerable 
importance  in  others,  for  example  in  assessing  flamability  and  the  health  impacts  of  toxic  release, 
in  studies  of  biological  response  to  odor  [14,  28],  and  pollution  damage  to  plants  [29].  Because 
behavioral  reactions  are  governed  by  instantaneous  or  nearly  instantaneous  concentrations  of  odor 
[30]  and  such  probability  density  estimates  cannot  be  used  to  accurately  predict  where  in  space  an 
odor  is  instantaneously  above  threshold,  probability  density  plume  models  are  not  appropriate  for 
the  CPT  strategy  performance  analysis. 

This  article  presents  a  model  for  dispersion  in  a  turbulent  medium.  The  specific  purpose  of 
this  plume  dispersion  model  is  to  facilitate  modeling  and  analysis  of  navigation  strategies  designed 


to  locate  a  wind-dispersed  odor  plume  and  then  to  trace  it  to  its  source^.  It  is  well  known  (e.g., 
[31,  32])  that  direct  numeric  simulation  (DNS)  of  turbulent  flow  in  reasonable  amounts  of  time  is 
limited  by  computational  power  to  large  grid  size  or  low  Reynolds  numbers.  Therefore,  this  article 
presents  the  design  of  a  simplified  plume  simulation,  based  on  physical  principles,  that  resulted 
in  a  computationally  feasible  simulation  (i.e.,  100  batches  of  100  source  relative  starting  locations 
evaluated  on  a  100  by  100  m  area  in  less  that  24  hours  on  a  typical  desktop  computer).  The 
simulation  design  presented  herein  addresses  the  following  tradeoffs: 

1.  Field  studies  have  shown  that  odor  detection  events  can  have  pulse  widths  on  the  order  of  0.01s 
[27].  Moths  have  been  shown  to  respond  to  odorant  pulse  widths  on  the  order  of  0.1  s  [14]. 
Therefore,  the  fine-time  scale  structure  of  the  sensed  odor  is  important. 

2.  Although  the  fine  structure  of  the  turbulent  flow  is  important  in  the  process  of  transporting 
relatively  undiluted  parcels  of  odor  over  long  distances,  the  anticipated  vehicle  flow  sensors  will 
be  unable  to  resolve  the  fine-scale  flow  structure.  Therefore,  accurate  simulation  of  the  fine-scale 
flow  itself  is  not  important. 

3.  The  plume  must  have  an  intermittent  internal  structure  in  the  sense  that  a  snapshot  of  the 
plume  appears  as  patches  of  above  threshold  odor. 

4.  The  sinuous  (or  meandering)  nature  of  the  plume  is  a  key  complicating  factor  to  plume  trac¬ 
ing.  The  simulation  must  generate  plumes  that  meander.  In  addition,  the  meander  should  be 
coherent  with  the  flow  field  in  the  sense  that  the  odor  distribution  down  wind  from  the  source 
is  the  result  of  advection  by  the  flow  field. 

5.  The  simulation  must  generate  exemplar  plumes,  not  time-averaged  odor  distributions.  The 
simulation  should  be  capable  of  returning  measured  concentration  at  location  x{t)  at  time  t. 
Both  position  and  time  will  change  as  the  searcher  executes  its  mission.  The  searcher  trajectory 
will  respond  to  the  sensed  flow  and  odor  detection  events.  This  searcher  trajectory  cannot  be 
predicted  prior  to  the  simulation.  Therefore,  the  simulation  must  maintain  models  of  the  flow 
and  odor  concentration  as  a  function  of  position  and  time. 

6.  The  effect  of  the  searcher  on  the  flow  and  the  odor  distribution  are  not  considered  important 
for  this  application,  because  the  searcher,  which  is  trying  to  make  progress  up  the  plume. 
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mainly  affects  the  downwind  flow  and  odor  distribution.  Inclusion  of  these  effects  would  signifi¬ 
cantly  increase  the  computational  burden  without  significantly  affecting  the  strategy  evaluation 
results. 

The  goal  of  the  simulation  discussed  herein  is,  therefore,  to  produce  a  more  challenging  and  physi¬ 
cally  plausible  simulation  model  than  those  currently  used  in  insect  behavioral  and  robotic  strategy 
design  studies,  e.g.,  [34,  33,  35,  36],  while  achieving  significant  computational  simplification  relative 
to  the  turbulence  simulation  models  in  the  literature,  e.g.,  [37,  38,  39,  40,  41,  32]. 

1.1.  Article  Overview 

To  generate  realistic  simulated  plumes  we  considered  the  following  three  characteristic  structural 
features  to  be  of  primary  concern.  First,  the  sensed  plume  at  a  fixed  location  should  have  an 
intermittent  internal  structure  that  closely  duplicates  experimental  observations.  Second,  the  plume 
that  results  from  a  semi-continuous  release  of  odor  should  be  sinuous  and  time  varying.  Third,  the 
plume  shape  and  wind  field  history  should  be  coherent. 

The  first  characteristic  feature  is  achieved  herein  by  representing  the  odor  plume  as  a  sequence 
of  puffs  [33,  55].  Puffs  are  released  sequentially  at  the  source  location.  Each  puff  is  composed  of 
n  filaments.  After  release,  the  filament  center  location  is  determined  by  integration  of  the  wind 
velocity: 

Pi  =  v(pi)  (2) 

where  pi  and  v(pi)  are  the  3  dimensional  position  of  the  i-th  filament  and  wind  velocity  at  the 
location  of  the  i-th  filament.  To  achieve  the  second  characteristic  feature,  our  model  incorporates 
a  dynamic  model  of  the  wind.  This  model  allows  the  wind,  which  transports  the  filaments,  to  be  a 
continuous  but  varying  function  of  position  and  time.  Changes  in  wind  direction  cause  the  plume 
to  meander  (form  a  snake-like  path  when  viewed  from  above).  As  the  wind  velocity  changes,  the 
instantaneous  wind  direction  at  given  positions  within  the  plume’s  boundaries  sometimes  will  not 
point  either  toward  the  odor’s  source  nor  be  coincident  with  the  direction  of  the  plume’s  centerline. 
These  characteristics  of  the  simulated  plume  reproduce  experimentally  observed  phenomena  [42,  43]. 

The  puff-based  model  presented  here  is  general  and  is  designed  so  that  it  can  be  tuned  for  a  range 
of  applications.  For  the  simulation  results  presented  here,  the  parameters  were  tuned  to  correspond 


to  odor  plumes  as  perceived  by  male  moths.  The  receptor  system  of  some  male  moths  is  capable 
of  resolving  10  odor  pulses  per  second  [44].  A  flying  male  can  react  to  odor  filaments  of  200  ms 
duration  [14,  45].  Therefore,  we  designed  the  present  model  of  plume  dispersion  with  an  integration 
time  step  of  10  ms.  Similarly,  the  units  of  the  sensed  odor  (see  Section  2.4)  were  specified  to  be 
relevant  to  the  male  moth. 

Conventional  time-averaged  concentration  models  are  in  general  time  independent.  It  is  assumed 
that  the  plume  was  fully  developed  prior  to  the  initiation  of  data  collection.  Such  concentration 
distributions  (assuming  stationarity  and  ergodicity)  allow  time  independent  predictions  over  o,n 
ensemble  of  plumes  of  the  average  concentration  as  a  function  of  position.  Measured  instantaneous 
concentration  of  specific  exemplar  plumes  can,  however,  deviate  significantly  from  the  ensemble 
average  [27,  46,  47,  48].  Alternatively,  instantaneous  concentrations  as  produced  by  simulation 
models  produce  a  model  instance  of  an  odor  plume  as  an  individual  realization.  By  design,  such  a 
plume  will  be  time  varying;  however,  when  time  averaged,  the  mean  concentration  should  match 
the  above  ensemble  averages.  In  addition,  simulated  plume  models  enable  certain  types  of  analysis 
that  are  precluded  by  the  loss  of  time  signature  information  in  mean  concentration  plume  models, 
such  as  analysis  of  alternative  biological  strategies  under  uniform  experimental  conditions  (see  e.g., 
[34,  33]). 

The  objective  of  this  article  is  to  present  a  simulation  model  for  an  odor  plume  that  (1)  realisti¬ 
cally  reproduces  the  short  time-scale  signature  for  sensed  concentrations,  (2)  accurately  reproduces 
the  long-term  time-averaged  plume  data,  (3)  incorporates  a  simulated  wind  (advection)  model  that 
is  continuous  over  a  region  of  interest  but  time  varying  in  magnitude  and  direction,  and  (4)  is 
computable  in  reasonable  time.  The  presentation  of  the  model  is  followed  by  analysis  of  simulated 
results  in  relation  to  experimental  results  previously  reported  in  the  literature. 

2.  Model  Overview 

A  chemical  released  at  a  source  location  will  be  manipulated  by  turbulent  and  molecular  diffusion 
while  transported  advectively  by  the  wind.  As  described  in  [46,  49],  the  odor  dispersion  process  is 
dominated  by  turbulent  dispersion.  Turbulent  dispersion  involves  a  wide  range  of  length  (or  eddy) 
scales.  Eddies  larger  than  the  scale  of  a  section  of  the  plume  (i.e.,  a  puff)  transport  the  puff  as  a 


whole,  causing  the  ensemble  of  puffs  to  appear  as  a  sinuous  plume.  Eddies  smaller  that  the  puff 
mix  the  components  of  the  puff  causing  little  puff  motion  or  growth.  Eddies  on  the  order  of  the 
puff  size  cause  significant  growth/distortion  of  the  puff  and  motion  of  individual  filaments  relative 
to  the  instantaneous  plume  centerline. 

To  develop  a  corresponding  plume  simulation  model,  the  velocity  vector  will  be  decomposed  into 
three  components:  v^,  \rn,  and  v^.  Each  component  will  be  implemented  by  a  distinct  process  as 
described  in  the  subsequent  sections.  This  decomposition  of  the  velocity  spectrum  is  motivated  by 
the  numerical  implementation,  but  can  be  interpreted  theoretically  in  terms  of  the  eddy  scales  of 
the  previous  paragraph.  To  clarify  this  correspondence,  it  is  first  necessary  to  state  that  the  puffs, 
for  simulation,  will  be  further  decomposed  into  parcels  referred  to  herein  as  filaments.  The  effect  of 
the  smallest  eddies  (i.e.,  slow  growth  of  the  filaments)  of  the  wind  fluid  flow  process  (modeled  by  v^) 
will  be  implemented  as  an  increase  in  filament  size  and  change  in  shape,  as  described  in  Section  2.3. 
The  term  represents  the  portion  of  the  wind  process  with  characteristic  length  much  larger  than 
the  filaments.  This  portion  of  the  wind  process  transports  each  filament  as  a  body;  therefore,  the 
term  v^,  represents  advection.  The  advection  portion  of  the  velocity  is  represented  as  a  continuous 
(in  time  and  space),  but  temporally  and  spatially  varying  function  (see  Section  2.1),  so  that  a 
sequence  of  filaments  released  at  the  source  will  result  in  a  sinuous  trail  of  filaments  leaving  the 
source.  The  term  represents  the  intermediate  range  of  scales  that  transports  (i.e.,  stirs)  the 
filaments  within  the  body  of  the  plume. 

Although  the  position  of  individual  molecules  cannot  be  conveniently  simulated,  the  idea  is 
useful  to  consider  for  the  purpose  of  defining  the  simulation.  Therefore,  the  position  of  a  single 
molecule  is  given  by 


Pm  =  Vrf  +  Vm  +  Va-  (3) 

Because  the  intermediate  and  advective  terms  affect  all  molecules  within  a  filament  similarly,  the 
model  can  be  decomposed  into  two  processes-the  changing  shape  of  the  filament  and  the  transport 
of  the  filament.  The  position  of  the  i-th  filament  will  be  represented  by  pj.  The  changing  position 
of  the  i-th  filament  is  represented  by 


Pi  =  Vmj  -I-  Va. 


(4) 


where  the  implementation  of  the  Vm;  and  Va  processes  is  discussed  in  subsequent  sections. 


Physically,  the  growing  and  changing  shape  of  the  filament  is  a  very  interesting  process  dom¬ 
inated  by  turbulent  diffusion.  However,  numerical  simulation  of  diffusion  processes  is  a  complex 
and  (simulation)  time  consuming  task.  Therefore,  at  this  level  some  compromise  is  required  to 
achieve  reasonable  simulation  times  while  still  achieving  the  objectives  of  the  simulation  model. 
The  approach  herein  is  to  use  a  template  filament  shape  that  can  change  size  and  orientation.  The 
simulation  results  are  not  sensitive  to  the  choice  of  template  shape,  as  long  as  the  change  in  shape 
is  physically  reasonable. 

Any  numeric  flow  simulation  must  satisfy  various  conservation  equations  (see  Ch.  2  in  [39]). 
For  the  simulation  herein,  conservation  of  the  mass  of  the  flowing  material  and  conservation  of 
momentum  is  addressed  in  Sections  2.1  and  2.2.  Conservation  of  the  mass  of  the  odor  is  addressed 
through  the  filament  representation  as  described  in  Section  2.3. 


2.1.  Advection 


This  section  describes  the  theoretical  model  for  =  {u,v,w)  suitable  for  simulation  implemen¬ 
tation.  The  advection  model  is  currently  implemented  in  two  dimensions  (i.e.,  x  and  y)  with  the 
assumption  that  for  the  area  of  interest  is  independent  of  z.  This  assumption  is  motivated  by  the 
assumption  that  the  plume  tracing  will  occur  at  a  nearly  constant  altitude  within  a  few  meters  of 
the  ocean  floor  or  earth  surface.  Issues  related  to  plume  dispersion  in  a  stratified  flow  are  discussed 
in  [50]. 

Within  the  atmospheric  surface  layer,  we  have  assumed  that  compared  with  turbulent  forces,  the 
following  forces  on  an  air  parcel  are  negligible:  Coriolis,  geostrophic  winds,  and  molecular  viscosity. 
It  is  also  assumed  that  time-averaged  pressure  terms  can  be  neglected,  so  that  the  time-mean 
motion  equations  are 


du 

_du 

_du 

_  du 

dt 

dz 

du^v! 

du'v^ 

du^w^ 

dx 

dy 

dz 

dv 

_  dv 

_dv 

_  dv 

dt 

"^Vz 

dv^u^ 

dv’v^ 

dv^w^ 

dx 

dy 

dz 

(5) 


(6) 
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with  the  continuity  equation 


du  dv  dw 
dx^  dz 


(7) 


where  u{t)  -  u{t)  -I-  u'{t).  To  conserve  computational  resources  we  have  used  the  simplest  form  of 
K-closure  method  [51,  39]  under  the  assumption,  applicable  under  near  neutral  conditions,  that 


u'u'  =  —Kx 


u'w'  =  -K, 


du 

dx 

du 

dz 


v'v'  =  -Ky^ 
dv 


and 


u'v>  =  v>u'  =  -l(^K.^+Ky^ 


(8) 

(9) 

(10) 


where  the  Kx,Ky,  and  Kz  terms  represent  difFusivity.  Assuming  that,  for  the  region  of  the  plume 
with  which  we  are  concerned  (approximately  a  hundred  meters  in  length),  turbulence  is  homoge¬ 
neous  and  approximately  isotropic. 


ATx  -  Ky  (11) 

=  ^  =  (12) 

dx  dy  dz 

Note  that  the  Kx.K^j,  and  Kz  terms  could  be  assumed  to  be  appropriate  functions  of  position  at 
minimal  computational  cost,  with  appropriate  changes  in  eqns.  (11-12).  In  the  present  simulation, 
Kx  G  [1,30]^  is  a  constant.  Under  these  assumptions,  the  dynamic  wind  (advection)  model  for 
the  u-term  simplifies  as  follows 


du 

dt 
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du 
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The  model  for  v  can  be  processed  similarly.  If  it  is  assumed  that  the  model  is  applied  in  a  narrow 
layer  near  the  surface  of  the  earth,  on  a  region  where  =  [u,v,w)  is  independent  of  then  the 
following  simplified  two  dimensional  model  results: 

du 


dv 
dt 

These  equations  can  be  solved  numerically  at  a  grid  of  points  on  the  region  of  interest  for  a  specific 
set  of  time-varying  boundary  conditions.  The  term  Vq  derived  in  this  section  has  the  sole  function 
of  transporting  the  fine-scale  (filaments)  of  the  plume.  Therefore,  the  grid  point  separation  A 
is  defined  to  be  larger  than  the  scale  of  the  filaments.  The  advective  wind  velocity  at  locations 
between  the  grid  points  is  determined  by  interpolation  based  on  the  advective  wind  at  the  adjacent 
grid  points.  This  approach  yields  a  continuous,  spatially  varying,  and  time  varying  model  of  the 
advection  of  the  plume  elements  with  a  reasonable  characterization  of  the  physics.  Therefore,  the 
wind  velocities  at  distinct  but  nearby  points  are  correlated,  but  spatially  varying. 

When  a  sequence  of  filaments  is  released  at  a  fixed  location,  the  advection  terms  (u,  v,  w)  generate 
a  meandering  trail  of  filaments,  as  shown  in  Fig.  1.  This  figure  depicts  a  top  down  view  of  a  100  by 
100  m  field.  The  above  equations  were  solved  numerically  for  boundary  conditions  generated  by  a 
mean  flow  of  [1,0]—  plus  a  colored  noise  process.  The  colored  noise  process  at  each  corner  is  im¬ 
plemented  by  filtering  unit  Gaussian  white  noise  by  H{s)  =  For  simulations  in  which  we 

want  to  generate  large  amplitude  meander,  we  typically  select  (a,  6,  G)  =  (0.04^^,  0.04^,  20).  For 
simulations  in  which  we  want  to  generate  small  amplitude  meander,  we  typically  select  (a,  6,  G)  — 
(0.5^^,  0.1^,  3).  The  boundary  condition  along  the  edge  nodes  is  then  generated  by  interpolating 
between  the  values  at  the  two  adjacent  corners.  The  components  of  the  model  that  describe  diffusion 
of  the  plume  elements  relative  to  the  plume  centerline  and  the  growth  of  the  filaments  were  turned 
off  for  the  simulation  that  generated  Fig.  1,  so  that  the  figure  only  shows  an  instantaneous  example 
of  plume  centerline  meander.  Centerline  relative  filament  diffusion  and  filament  growth  will  be 
discussed  in  subsequent  sections.  The  arrows  on  the  figure  indicate  the  instantaneous  wind  vector 
at  the  tail  of  the  arrow.  The  instantaneous  wind  vector  and  the  local  plume  centerline  are  highly 


_du  _du 
1  1 
_dv  _dv 
1  d'^v  1  ^  d'^v 


(17) 


(18) 


correlated  near  the  source  of  the  plume.  As  the  distance  from  the  source  of  the  plume  increase,  the 
correlation  between  the  directions  of  these  two  vectors  decreases.  The  physical  explanation  of  such 
behavior  is  that  each  filament  is  transported  instantaneously  along  the  local  wind  vector;  however, 
the  plume  centerline  alignment  is  determined  by  the  entire  past  history  of  local  wind  experienced 
by  all  filaments  making  up  the  plume. 

Implementation  of  the  advective  model  as  described  herein  yields  a  spatially  and  time  varying 
flow  field  over  an  area  of  interest.  Transport  of  odor  filaments  by  this  flow  field  yields  a  plume 
centerline  with  realistic  characteristics.  In  particular,  the  meandering  nature  of  the  plume  (with 
‘semi-independent’  instantaneous  wind  direction  and  plume  centerline  alignment)  is  a  realistic 
characteristic  of  natural  plumes  [42,  43]  and  therefore  is  a  significant  advancement  relative  to 
the  wind-tunnel  type  of  puff  models  used  previously  for  analysis  of  insect  flight  response,  e.g.  see 
[34,  33].  The  finite  element  model  for  Va  is  not  intended  to  model  the  small-scale  turbulent  structure 
of  the  wind.  Therefore,  a  relatively  large  spatial  separation  (i.e.,  5-10  m)  of  the  finite  element  nodes 
is  reasonable.  This  large  separation  allows  the  approach  to  be  computationally  feasible. 


2.2.  Centerline  Relative  Filament  Diffusion 

This  section  describes  the  model  for  Vm-  This  is  the  scale  of  the  flow  between  A  and  the  filament 
size.  Assuming  that  A  is  specified  judiciously,  this  scale  of  the  flow  distributes  the  filaments  about 
the  plume  centerline.  For  computational  reasons  (i.e.,  the  number  of  nodes  increases  as  ^  for 
three  dimensional  implementations),  this  scale  of  the  flow  is  not  modeled  through  finite  difference 
methods.  Instead,  the  velocity  of  the  i-th  filament  relative  to  the  centerline  is  modeled  as  a 
random  process  implemented  in  state-space  notation  [52]  as 

U  =  A  (jO  B  U  ^rni  —  C  LO  D  V  (19) 

where  u)  €  A,  B,  C,  and  D  are  appropriately  sized  matrices;  and  z/  is  a  white  noise  process  with 
power  spectral  density  al  (i.e.,  cov  A  r))  =  alSir)  where  5(t)  is  the  impulse  function,  see 

p.  120  in  [52]  or  pp.  81-85  in  [53]).  The  transfer  function  from  i/  to  v^^,.  is 

^  =  H{s)  =  C{s  I-A)-'^B  +  D 

v 
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(20) 


where  s  is  the  Laplace  variable.  Therefore,  is  a  colored  noise  process  with  spectral  density 
defined  by 

PSD„^{s)  =  H*[s)H(s)al  (21) 

where  is  the  complex  conjugate  transpose  of  H. 

Note  that  if  the  bandwidth  of  H  is  large,  then  'Vrn  ^  where  G  =  if(0),  so  that  Vyyj  is  a 
white  noise  process.  In  this  case,  let  Y  denote  the  standard  deviation  of  the  filament  position  in 
the  y-direction  (i.e.,  plume  width)  due  only  to  (i.e.,  relative  to  the  centerline)  .  In  this  case, 
Y  —  G(Jtyy/f,,  where  t  is  the  time  since  the  puff  was  released.  This  matches  Robert’s  model  (see  [54] 
and  eqn.  (21)  in  [55]). 

Fig.  2  displays  an  example  of  the  plume  generated  by  the  relative  diffusion  terms  when 
advection  is  the  constant  vector  ==  [1,0,0]  ra/s  and  filament  growth  is  turned  off;  therefore,  this 
figure  illustrates  only  the  effect  of  Vyy^  (i.e.,  centerline  relative  diffusion).  For  this  simulation,  each 
component  of  was  a  white  noise  process  with  spectral  density  given  by  ai^  ==  2;^^. 

2.3.  Filament  Model  and  Growth 


When  summed,  and  v„.  describe  the  velocity  of  the  filament  center  as  the  filament  is  bodily 
transported  over  the  field  of  interest.  The  pointwise  calculation  of  concentration  and  changing  shape 
of  the  filament  are  the  topics  of  this  section.  As  shown  in  Fig.  2,  the  plume  is  the  composition 
of  a  large  number  of  advected  and  dispersed  filaments.  Given  a  large  number  of  filaments,  the 
overall  instantaneous  concentration  at  x  =  {x^y^z)  due  to  all  the  filaments  is  just  the  sum  of  the 
concentrations  at  that  location  contributed  by  each  filament: 


N 


2  =  1 


molecules 


cm'^ 


(22) 


where  N  is  the  number  of  filaments  currently  being  simulated.  Note  that  this  model  allows  the  in¬ 
stantaneous  odor  concentration  to  be  evaluated  at  any  position  and  time  of  interest.  This  pointwise 
concentration  measurement  can  be  converted  into  a  flux  by  multiplying  by  the  local  wind  velocity 
Vft  and  the  effective  area  of  the  sensor  cross-section. 

The  concentration  at  location  x  due  to  the  i-th  filament  is  modeled  as 

r  (  --  Q  ( molecules 

^(x,  exp  ^  ^2^^^  j  cm^ filament 

ri{t)  =  ||x~py(t)||  cm 
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where  Q  represents  the  amount  of  odor  released  (i.e.,  molecules  per  filament)  and  is  a  parameter 
controlling  the  size  of  the  i-th  filament.  This  expression  is  derived  based  on  the  assumption  that  the 
filament  contains  material  normally  distributed  relative  to  the  filament  center.  Therefore,  if  (7?:(x) 
is  integrated  over  the  (infinite)  spatial  extent  of  the  filament,  then  it  correctly  predicts  Q  molecules 
in  the  filament.  Therefore,  the  mass  of  the  odor  producing  chemical  is  conserved  as  the  filament 
size  Ri  changes. 

Various  models  are  available  for  the  time  derivative  of  the  radius  of  the  i-th  filament.  Here,  only 
two  models  are  presented.  For  the  first  model,  let  the  radius  of  the  i-th  filament  be  assumed  to 
change  as 

R(t)  =  +7^ 


then  the  rate  of  growth  of  the  i-th  filament  is 


(m 

dt 


Note  that  the  same  solution  can  be  obtained  with  fewer  computations  by  solving 


dp 


and  defining  R  =  p^'^.  For  the  second  model,  let  the  radius  of  the  i-th  filament  be  assumed  to 
change  as 

then  the  rate  of  growth  of  the  i-th  filament  is 

dR  _  j 
dt  2R 


Note  that  the  same  solution  can  be  obtained  with  fewer 

dp 


dt 


=  7 


computations  by  solving 


and  defining  R-p^.  These  expressions  are  intended  to  account  for  molecular  diffusion  and  growth 
of  the  filament  due  to  the  smallest  length  scales  of  the  turbulent  flow.  Alternative  models  of  filament 
growth  could  be  combined  to  account  for  different  regimes  of  filament  growth. 

For  the  simulations  presented  herein,  Q  =  8,3  x  10^  rnoiecuies  _  ^  pheromone  with 

molecular  weight  282  which  is  the  rate  of  emission  of  pheromone  by  a  female  gypsy  moth 

[56].  The  parameter  Q  —  ^  where  n  is  the  filament  release  rate  in  ,  The  filament  release 
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rate  is  of  considerable  importance  and  is  a  user  determined  function.  All  simulations  presented 
herein  have  n  =  10 The  second  model  of  filament  growth  is  used  with  7  =  0.001^  and 

=  O.OOlm^  at  the  time  of  release.  This  initial  filament  size  was  selected  to  be  characteristic  of 
the  size  of  the  female  gypsy  moth. 

Fig.  3  displays  an  example  of  the  plume  that  results  when  advection,  centerline  relative  diffu¬ 
sion,  and  filament  growth  are  all  active.  Again,  the  local  instantaneous  wind  direction  and  plume 
centerline  do  not  necessarily  align,  except  near  the  source.  Fig,  4  displays  a  typical  example  of  the 
sensed  concentration  resulting  from  the  simulated  model. 

An  alternative  expression 


C't(x)  = 


Q 


\/8n^  det(P) 


exp(-(x-pi)P  ^(x-pi)) 


molecules 


cm'^ 


(23) 


with  P  a  time  varying  positive  definite  matrix  would  allow  both  the  size  and  shape  of  each  filament 
to  be  altered  while  maintaining  the  number  of  molecules  per  filament  as  Q  (i.e.,  conserving  the 
mass  of  the  the  odor-producing  chemical). 


2.4.  Sensor  Model 


The  sensed  concentration  is  modeled  as  a  low  pass  filtered  and  threshold  specified  version  of  the 
instantaneous  concentration, 


dc{t) 

dt 

y{t)  = 


ac{t)  +  aC(xs(t)) 

[  c(t)  if  c{t)  >  T 
I  0  otherwise 


(24) 

(25) 


where  a  is  the  filter  bandwidth,  t  is  the  sensor  threshold,  c{t)  is  an  internal  state  of  the  filter, 
and  Xs(t)  is  the  (possibly  time  varying)  sensor  position.  The  filter  input  is  the  instantaneous 
concentration  at  the  sensor  location  C'(xs(t)).  The  filter  output  is  y{t).  The  effect  of  the  sensor 
on  the  downwind  plume  is  not  modeled.  In  the  results  that  follow,  a  is  varied.  The  threshold  r 
was  selected  to  match  the  threshold  of  male  response  to  pheromone  in  the  gypsy  moth  {Lymantria 
dispar)  which  is  4  X  [30]. 
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3.  Model  Analysis 


This  section  presents  simulated  plume  data  in  various  formats  (long-term  averages,  amplitude 
statistics,  temporal  statistics).  The  intent  is  to  present  sufficient  information  to  allow  verification 
that  the  simulated  plume  has  the  general  characteristics  of  odor  plumes  encountered  in  field  ex¬ 
periments.  Our  main  interests  in  such  a  validation  are  the  amplitude  (i.e.,  mean  versus  distance, 
’in-plume’  mean  versus  distance,  and  peak  to  mean)  and  short  duration  temporal  statistics  (i.e., 
intermittency,  pulse  width,  and  burst  return  interval).  The  higher  order  statistical  information  is 
presented  for  completeness.  See  the  Appendix  for  definitions  of  the  various  statistical  quantities. 

For  analysis  relative  to  field  experiments,  we  compare  the  simulation  results  with  field  mea¬ 
surements  presented  by  Jones  in  [27].  Jones  analyzed  experimental  data  acquired  using  specially 
designed  instruments  with  a  time  resolution  of  0.01  s.  The  data  presented  by  Jones  includes  plume 
statistics  as  a  function  of  sensor  bandwidth  and  distance  from  the  source.  This  portion  of  the  Jones 
data  is  included  herein  for  ease  of  reference  as  Tables  I  and  11.  In  addition,  the  Jones  article  contains 
histograms  that  present  short-term  odor  statistics  (pulse  amplitude,  burst  length,  and  burst  return 
period)  that  are  directly  relevant  to  the  studies  motivating  the  development  of  this  simulation 
model.  Those  figures  could  not  be  included  herein,  but  they  are  compared  qualitatively  with  the 
simulation  results  in  Section  3.3. 

3.1.  Long- Duration  Time  Averages 

Gaussian  models  are  widely  used  descriptions  for  the  long-time  average  concentration.  Although 
our  main  interest  is  to  replicate  in  simulation  the  short-term  intermittent  structure  of  plumes,  it  is 
also  of  interest  to  investigate  whether  the  long-term  time  average  of  the  simulated  plume  replicates 
existing  Gaussian  models  [24]. 

Fig.  5  plots  the  3  minute  time-average  contour  for  the  simulated  plume  corresponding  to  the 
threshold  of  0.04  x  shown  is  the  Sutton  model  isopleth,  generated  from  eqn.  (1) 

with  Sy  =  O.SCyX^"^/^  and  Sz  —  for  the  parameters  n=l,  Q==20,  Cy=0A,  and  Cz=0.2. 

The  3-minute  time-averaged  simulation  data  contour  is  in  close  agreement  with  the  Gaussian 
contour.  This  is  true  even  though  the  instantaneous  plume  has  the  desired  characteristics  (i.e., 
sinuous  and  intermittent).  This  fact  is  highlighted  by  a  direct  comparison  of  Figs.  3  and  5.  In 


simulated  meters,  both  figures  have  the  same  x  dimensions.  The  y-scale  of  Fig.  3  is  twice  the  scale 
of  Fig.  5.  Fig.  3  plots  instantaneous  concentration  as  a  function  of  position,  whereas  Fig.  5  shows 
the  time-average  contour.  Misuse  of  Fig.  5  would  imply  that  a  sensor  placed  anywhere  within  the 
contour  would  detect  above  threshold  concentrations.  The  instantaneous  plume  ‘snapshot’  of  Figs. 
3  shows  that  the  above  threshold  locations  are  actually  localized  and  time  varying. 

For  clarity  of  presentation,  additional  time-averaged  simulated  plume  contours  are  not  included 
in  Fig.  5;  however,  such  contours  match  additional  characteristics  of  the  Gaussian  plume  model. 
For  example,  as  the  duration  of  the  time-average  increases,  the  width  of  a  given  contour  increases. 

3.2.  Amplitude  Statistics 

Statistics^  related  to  the  amplitude  of  the  concentration  time  series  are  presented  in  this  section. 
Because  the  intent  of  this  section  is  to  validate  that  the  simulation  data  reproduce  key  characteristics 
of  actual  plume  data,  the  form  of  the  data  presentation  herein  is  defined  to  match  that  of  [27].  To 
facilitate  comparison,  the  tabular  data  from  [27]  is  included  here  as  Tables  I  and  II.  The  simulated 
data  used  to  compute  the  statistics  of  this  section  are  from  a  simulation  with  an  integration  time 
step  of  0.01  s  and  a  duration  of  10  minutes  (i.e.,  60000  points).  During  the  simulation,  the  filament 
release  rate  was  10  filaments  per  second.  The  simulated  data  statistics  are  presented  in  Tables  III 
and  IV. 

The  first  row  of  data  in  Table  III  presents  mean  sensed  concentration  as  a  function  of  downwind 
distance  {d  =2,  5,  10,  and  15  meters).  As  for  the  Jones  data  in  Table  I,  the  mean  concentration 
decreases  as  a  function  of  distance  from  the  source.  The  top  graph  of  Figure  6  plots  the  mean 
concentration  versus  downwind  distance  from  the  source  for  the  simulation  and  Jones  data.  Because 
the  units  are  distinct,  both  sets  of  data  have  been  normalized  to  unit  magnitude  at  a  downwind 
distance  of  2  m.  The  decrease  of  the  mean  with  downwind  distance  is  caused  by  the  width  of  the 
plume  increasing,  each  puff  becoming  wider,  and  the  area  over  which  the  plume  meanders  increasing 
with  distance  from  the  source.  As  noted  in  [27],  the  mean  concentration  is  not  affected  by  filtering 
of  the  data. 

The  second  row  of  data  in  Table  III  presents  the  conditional  mean  of  the  simulated  plume  as 
a  function  of  downwind  distance  from  the  source.  This  quantity  is  determined  by  calculating  the 
mean  concentration  only  when  the  sensor  is  in  the  plume.  Note  that,  as  expected  based  on  physical 
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principles,  the  conditional  mean  is  larger  than  the  unconditional  mean  and  the  conditional  mean 
decreases  more  slowly  than  the  unconditional  mean  as  a  function  of  source  distance.  The  bottom 
graph  of  Figure  6  plots  the  conditional  ‘in  the  plume’  mean  concentration  versus  downwind  distance 
from  the  source  for  the  simulation  and  Jones  data.  Again  both  sets  of  data  have  been  normalized 
to  unit  magnitude  at  a  downwind  distance  of  2  m. 

Table  IV  presents  second  through  fourth  order  (non-dimensional)  moments  and  other  key  pa¬ 
rameters  of  the  concentration  time  series  as  a  function  of  both  downwind  distance  and  filtering 
bandwidth.  The  indicated  bandwidth  is  that  of  a  first  order  low  pass  filter  (20  dB  per  decade). 
The  data  in  Table  IV  is  presented  in  the  same  format  as  the  Jones  data  of  Table  11.  Although 
it  is  very  difficult  to  match  higher  order  statistics,  the  general  trends  of  the  data  in  these  tables 
are  very  similar.  The  peak-to-mean  ratio  decreases  as  the  filter  bandwidth  decreases,  because  the 
peak  of  the  filtered  quantity  is  decreased  while  the  mean  is  unaffected.  Also,  the  peak-to-mean 
ratio  increases  with  distance,  because  the  mean  concentration  decreases  more  rapidly  (decreased 
by  filament  growth,  dispersion  about  the  centerline,  and  movement  of  the  centerline)  than  does  the 
peak  concentration  (decreased  by  filament  growth).  Lastly,  the  intermittency  (percent  out  of  the 
plume)  increases  with  distance,  because  the  meandering  of  the  plume  causes  the  plume  to  move 
across  a  larger  area  at  greater  distances  from  the  'source. 

3.3.  Temporal  Statistics 

Temporal  statistics  related  to  the  concentration  time  series  are  analyzed  and  presented  in  this 
section.  Again,  comparison  is  made  with  the  experimental  statistics  presented  in  [27].  The  main 
temporal  statistics  that  are  of  interest  are  the  distribution  of  burst  lengths  and  the  distribution  of 
burst  return  periods,  where  each  of  these  terms  is  defined  in  the  appendix.  These  distributions  are 
presented  in  Fig.  7  for  two  downwind  distances  and  for  two  thresholds.  The  threshold  is  important, 
as  it  defines  when  a  sensed  concentration  is  counted  for  constructing  the  histogram.  The  dynamic 
range  of  the  sensor  and  data  in  this  experiment  was  Cs{t)  G  [0.04, 10000]  x  10^  ^ 

The  histograms  correspond  well  with  those  of  [27].  In  both  the  field  and  simulation  histograms, 
the  burst  lengths  are  significantly  shorter  than  1.0  s  as  dictated  by  the  filamentous  nature  of  the  real 
and  simulated  plumes.  As  the  threshold  increases,  the  number  and  duration  of  the  bursts  decreases. 
In  both  the  field  and  simulation  data,  the  burst  return  distribution  appears  to  be  bimodal  with 
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one  peak  at  or  below  0.1  s  and  a  second  peak  at  or  above  1.0  s.  The  physical  explanation  for  this 
is  that  the  sensor  may  not  detect  concentrations  either  due  to  the  sensor  being  out  of  the  plume  or 
the  sensor  being  between  filaments  while  still  ‘in  the  plume.’  Meandering  of  the  plume  causes  the 
‘out  of  plume’  burst  return  period  to  have  long  duration  and  to  increase  with  distance  from  the 
source.  The  burst  return  period  between  filaments  while  ‘in  the  plume’  has  much  shorter  duration 
and  increases  slowly  with  distance. 


4.  Conclusions 

This  article  has  presented  a  computationally  feasible  simulation  model  suitable  for  Monte  Carlo  type 
analysis  of  dispersion  in  a  turbulent  medium.  This  tool  is  particularly  important  in  studies  where 
the  critical  item  of  interest  is  instantaneous  concentration  rather  than  time-averaged  exposure. 
The  method  presented  paid  particular  attention  to  creating  a  model  with  a  continuous  but  time 
varying  plume  centerline,  dispersion  about  the  centerline,  experimentally  valid  long-term  averages, 
and  experimentally  validated  short  duration  temporal  and  amplitude  statistics.  The  model  allows 
calculation  of  the  wind  velocity  and  concentration  at  any  location  and  at  any  simulated  time  step 
(increments  of  0.01  s  for  the  current  implementation).  The  user  definable  parameters  allow  the 
model  to  be  tuned  for  a  wide  variety  of  applications. 

To  achieve  the  fast  numeric  computation  required  for  the  batch  Monte  Carlo  simulation  analysis 
of  possible  plume  tracing  strategies,  the  plume  simulation  did  not  attempt  to  model  the  fine-scale 
fiow  characteristics.  This  simplifying  assumption  was  motivated  by  the  characteristics  of  the  plume 
tracing  application,  where  the  vehicle  is  expected  to  have  only  a  low  bandwidth  («1  Hz)  flow  sensor. 
The  simulation  was  designed  to  maintain  the  plume  characteristics  that  significantly  complicate 
the  plume  tracing  problem  (meander,  intermittency,  and  varying  flow).  The  simulation  therefore 
provides  a  challenging  environment  for  evaluation  of  plume  tracing  algorithms. 

The  level  of  fine-scale  resolution  that  is  used  by  biological  entities  or  that  is  required  for  engi¬ 
neered  devices  to  perform  plume  tracing  is  still  an  open  issue^.  For  biological  entities  the  answer 
to  this  question  is  a  key  issue  in  research  focused  on  identifying  how  nature  solves  such  problems. 
For  engineered  systems,  the  answer  to  this  question  determines,  in  part,  the  quality  of  sensors 
that  will  be  required.  If  the  task  can  be  solved  without  detailed  analysis  of  the  fine-scale  plume 
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structure,  then  low  performance  and  hence  lower  cost  sensors  can  be  used.  If  detailed  analysis  of 
the  fine  scale  structure  is  required,  then  simulation  based  performance  analysis  will  require  higher 
fidelity  modeling.  Some  possible  directions  for  improving  the  simulation  include;  (1)  adapting  the 
sub-gridscale  velocity  characteristics  based  on  the  characteristics  of  the  advective  velocity  [57,  58]; 
and,  (2)  adjusting  the  characteristics  of  v-mi  as  a  function  of  the  scale  of  the  puff  (i.e.,  Richardson  or 
scale-dependent  dispersion).  A  key  issue  related  to  such  improvements  will  be  the  tradeoff  between 
increased  model  fidelity  versus  increased  computation. 
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Appendix 

A.  Parameter  Definitions 

Tables  II  and  IV  presents  various  statistics  of  the  simulated  plume.  This  appendix  defines  the 
presented  statistical  quantities. 

A.l.  Amplitude  Statistics 

The  concentration  C(p,i)  is  a  time  varying  function  of  position.  Although  the  concentration  is 
time  varying,  we  will  assume  that  the  statistics  of  the  process  are  stationary  and  ergotic.  Therefore, 
ensemble  statistics  will  be  calculated  based  through  time  integrals. 

The  mean  concentration  is  defined  by 

f(p)  =  ^^  Cip,T)dT.  (26) 
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The  n-th  central  moment  is  defined  as 


M, 


.(P)  =  ^^  [C'(P,'^) -r(p)]"dr. 


(27) 


For  n=2,3,4,  these  moments  are  the  fluctuation  variance,  skev/ness,  and  kurtosis,  respectively.  To 
allow  direct  comparisons  between  different  experiments,  it  is  useful  to  express  the  n-th  central 
moment  in  the  non-dimensional  form 


-W„(p)  = 


(i/o^[C(p.T)-r(p)]”dr)'^’ 

r(p) 


(28) 


The  non-dimensional  second  moment  is  referred  to  as  the  relative  intensity  of  fluctuations.  In  Table 
IV  the  skewness  and  kurtosis  are  expressed  in  non-dimensional  form. 

The  quantity  f/f  is  the  peak  to  mean  ratio.  This  quantity  is  calculated  herein  as  the  peak  over 
the  entire  experiment  divided  by  the  mean  over  the  entire  experiment.  This  quantity  could  also 
be  calculated  as  the  average  peak  height  (over  all  peaks)  divided  by  the  mean  concentration  or 
as  the  average  peak  height  (over  all  peaks)  divided  by  the  conditional  (i.e.,  in  the  plume)  mean 
concentration. 

Due  to  the  sinuous  nature  of  the  plume,  there  exist  periods  of  time  when  the  sensor  will  not  be 
‘in  the  plume.’  Relative  to  a  threshold  t,  the  intermittency  I  is  defined  as  the  percentage  of  time 
during  the  experiment  that  the  sensed  concentration  was  below  threshold.  The  above  amplitude 
statistics  can  be  calculated  over  the  entire  duration  of  the  experiment  or  only  during  the  portion 
of  the  experiment  when  the  sensor  was  in  the  plume.  In  this  article,  only  the  second  row  of  data  in 
Tables  I  and  III  is  restricted  to  times  when  the  sensor  was  in  the  plume. 


A.2.  Temporal  Statistics 

Fig.  4  displays  a  short  segment  of  the  simulated  concentration  at  a  position  2  m  downwind  from 
the  source.  The  symbol  tpi  defines  the  width  of  the  i-th  peak  (i.e.,  above  threshold  t).  The  symbol 
tgi  defines  the  width  of  the  i-th  gap  (i.e.,  below  threshold  r).  The  symbol  4,  =  tpi  4-  tgi  defines  the 
time  between  the  start  of  two  consecutive  pulses  (i.e.,  pulse  return).  Note  that  these  time  periods 
and  their  distribution  are  threshold  dependent.  Fig.  4  indicates  the  peak  and  gap  duration  of  the 
first  pulse  corresponding  to  a  threshold  of  2000.  Based  on  these  definitions,  if  ripi  and  rigi  denote 
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(29) 


the  number  of  pulses  and  gaps  in  the  experiment,  respectively,  then 

Tipi 

2=1  2=1 

where  T  is  the  experiment  duration  and  the  approximation  is  only  due  to  the  experiment  not 
starting  and  stopping  at  the  beginning  or  end  of  pulses.  Fig.  7  displays  normalized  histograms  of 
the  pulse  widths  (i.e.,  tpi)  and  pulse  returns  (i.e.,  4i), 


Notes 

^  In  the  last  decade  there  have  been  at  least  three  research  programs  in  the  United  States:  the  DARPA  Dog’s  Nose 
Program,  the  ONR  Chemical  Sensing  in  the  Marine  Environment  Program,  and  the  DARPA/ONR  Chemical  Plume 
Tracing  Program. 

^  A  typical  Monte  Carlo  batch  evaluation  of  a  single  strategy  involves  10,000  (100  repetitions  of  each  starting 
location  from  a  grid  of  100  source  relative  starting  locations.)  plume  simulations  each  300  seconds  long.  A  single 
plume  simulation  must  be  capable  of  returning  at  time  t  the  concentration  and  (low  bandwidth)  flow  velocity  at  any 
location  within  a  100  x  100  m  search  area. 

^  A  few  simulation  images  are  contained  in  this  article.  The  time  varying  simulation  output  is  more  interesting. 
Therefore,  an  executable  version  (for  Windows)  of  the  simulation  described  herein  is  available  for  download  at 
http://www.ee. ucr.edu/^farrell.  The  source  code  in  C  is  available  through  the  first  author. 

^  See  appendix  for  parameter  definitions. 

^  The  authors  also  thank  the  reviewers  for  their  time  and  efforts.  The  ideas  of  this  paragraph  were  formed  based 
on  suggestions  from  the  reviewers. 
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(a;,  y)  =  (5,  0)m.  Each  arrow  indicates  the  local  wind  vector  at  the  tail  of  the  arrow.  The  meandering  of  the  plume 
centerline  is  due  to  advection  of  a  sequence  of  odor  filaments  by  the  local  wind. 


:0.0,-50.0) 


(100.0.-50.0; 


represented  is  100  by  100  m  with  the  odor  source  at  {x,y)  =  (5,0)m,  Each  arrow  indicates  the  local  wind  vector  at 
the  tail  of  the  arrow.  The  centerline  of  the  instantaneous  plume  does  not  meander  because  the  wind  is  constant. 


Figure  3.  Meandering  plume  with  centerline  relative  diffusion  and  filament  growth.  The  area  represented  is  100  by 
100  m  with  the  odor  source  at  [x^y)  =  (5,0)m.  Each  arrow  indicates  the  local  wind  vector  at  the  tail  of  the  arrow. 


Concentration.  10°  Molecules/cm 


Figure  4.  Concentration  time  series  2  m  downwind  from  the  source.  The  time  resolution  is  0.01  sec.  Release  rate  is 
3  pufFs/sec.  The  concentration  is  expressed  in  10® 
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Figure  5.  Contour  plots  for  Sutton  model  and  a  three  minute  average  of  the  simulated  sensed  concentration 
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Distance  downwind,  m 

Figure  6,  Plots  comparing  the  mean  concentration  and  4n  plume’  mean  concentration  data  between  simulation  (solid 
with  ‘x’)  and  field  data  from  [27]  (solid  with  The  actual  data  contained  in  Tables  I  and  III  has  been  normalized 
by  their  values  at  2.0  m. 


2m  downwind 


15m  downwind 


Figure  7.  Top  left — Burst  length  distributions  for  a  stationary  sensor  located  2  m  downwind  from  the  source.  Top 
right — Burst  length  distributions  for  a  stationary  sensor  located  15  m  downwind  from  the  source.  Bottom  left  Burst 
return  period  distributions  for  a  stationary  sensor  located  2  m  downwind  from  the  source.  Bottom  right  Burst  return 
period  distributions  for  a  stationary  sensor  located  15  m  downwind  from  the  source.  Compare  with  Figs.  6  and  7  in 


Table  1.  Data  from  [27].  The  mean  over  all  data  is  presented  in  the  top 
data  line.  The  conditional  mean  for  Hn  the  plume  data’  is  presented  in 
the  second  data  line. 


Downwind  distance  (m) 

2 

5 

10 

15 

Mean  Concentration,  f  (^) 

4.21 

0.525 

0.239 

0.159 

Conditional  Mean  Concentration,  F 

28.4 

5.33 

1.46 

0.549 

Table  II.  Data  from  [27].  Plume  data  statistics  as  a  function  of  downwind  distance  and 
sensor  bandwidth.  Processed  data  for  2  m  downwind  is  presented  in  the  first  five  data  rows. 
Processed  data  for  5  m  downwind  is  presented  in  the  second  five  data  rows.  Processed 
data  for  10  m  downwind  is  presented  in  the  final  five  data  rows.  Skewness  and  kurtosis 
are  expressed  in  non-dimensional  form. 


With  low-  pass  filtering 


Name 

Symbol 

Unfiltered 

30  Hz 

10  Hz 

3  Hz 

1  Hz 

0.3  Hz 

Standard  deviation 

12.6 

11.0 

8.71 

5.60 

3.77 

2.40 

Skewness 

‘S'r 

4.95 

4.87 

4.62 

4.31 

4.18 

3.93 

Kurtosis 

Kr 

30.2 

29.6 

27.6 

25.5 

23.2 

21.5 

Peak  to  mean  ratio 

f/f 

36.4 

31.9 

25.6 

18.4 

13.9 

6.74 

Intermittency 

n%) 

85.2 

82.7 

80.9 

78.8 

79.1 

86.5 

Standard  deviation 

2.21 

2.08 

1.83 

1.39 

1.03 

0.722 

Skewness 

Sr 

7.18 

6.46 

6-27 

5.47 

4.24 

3.29 

Kurtosis 

Kr 

66.6 

53.3 

51.4 

39.3 

23.7 

13.0 

Peak  to  mean  ratio 

f/f 

78.2 

62.7 

58.9 

43.3 

22.2 

10.6 

Intermittency 

m 

90.1 

87.9 

82.6 

84.8 

81.0 

82.9 

Standard  deviation 

0.818 

0.765 

0.679 

0.536 

0.395 

0.279 

Skewness 

Sr 

8.82 

7.82 

7.48 

5.62 

4.49 

3.50 

Kurtosis 

Kr 

129 

97.0 

89.8 

47.0 

30.6 

14.9 

Peak  to  mean  ratio 

f/f 

112 

90.4 

79.0 

48.0 

28.5 

12.2 

Intermittency 

/(%) 

83.7 

83.9 

85.2 

83.7 

83.7 

85.2 

34 


Table  III.  Mean  concentration  of  simulation  data  over  10  minutes  (60000 
time  samples)  as  a  function  of  distance  from  the  source  (unfiltered  data). 
The  mean  over  all  data  is  presented  in  the  top  data  line.  The  conditional 
mean  for  ‘in  the  plume  data’  is  presented  in  the  second  data  line.  This  table 
is  compared  with  the  field  data  of  Table  I  in  Fig.  6. 


Downwind  distance  (m) 

2 

5 

10 

15 

Mean  Concentration,  f 

698.01 

232.22 

69.38 

40.48 

Conditional  Mean  Concentration,  f 

1494.0 

655.4 

237.3 

136.9 

Table  IV.  Simulated  plume  data  statistics  as  a  function  of  downwind  distance  and  sensor  band¬ 
width.  Processed  data  for  2  m  downwind  is  presented  in  the  first  six  data  rows.  Processed  data 
for  5  m  downwind  is  presented  in  the  second  six  data  rows.  Processed  data  for  10  m  downwind 
is  presented  in  the  final  six  data  rows.  Ten  minutes  (60000  samples)  of  data.  Compare  with  field 
data  of  [27]  that  is  contained  in  Table  II.  Skewness  and  kurtosis  are  expressed  in  non-dimensional 
form. 


With  low-  pass  filtering 


Name 

Symbol 

Unfiltered 

30  Hz 

10  Hz 

3  Hz 

1  Hz 

0.3  Hz 

Standard  deviation 

_  ji/finolec. 

2064 

2037 

1923 

1608 

1206 

873 

Skewness 

Sr 

4.53 

4.45 

4.14 

3.36 

2.32 

1.52 

Kurtosis 

Kv 

5.81 

5.70 

5.29 

4.21 

2.89 

1.80 

Peak  to  mean  ratio 

f/f 

14.33 

14.33 

14.33 

14.03 

11.76 

7.49 

Intermittency 

n%) 

53.28 

52.58 

48.54 

38.27 

28.11 

13.15 

Relative  intensity 

7r 

2.96 

2.92 

2.76 

2.30 

1.73 

1.25 

Standard  deviation 

<Tr. 

958 

944 

889 

750 

570 

411 

Skewness 

Sr 

7.75 

7.57 

6.90 

5.47 

3.80 

2.46 

Kurtosis 

Kr 

11.20 

10.88 

9.78 

7.54 

5.08 

3.15 

Peak  to  mean  ratio 

f/f 

43.06 

43.06 

42.59 

36.63 

23.58 

13.73 

Intermittency 

/(%) 

64.57 

64.38 

63.00 

56.84 

47.56 

28.70 

Relative  intensity 

7r 

4.13 

4.07 

3.83 

3-23 

2.46 

1.77 

Standard  deviation 

339 

333 

312 

265 

204 

148 

Skewness 

Sr 

11.13 

10.66 

9.40 

7.30 

5.02 

3.15 

Kurtosis 

Kr 

18.82 

17.76 

15.07 

11.13 

7.21 

4.17 

Peak  to  mean  ratio 

f/f 

144.14 

143.88 

130.83 

81.64 

46.86 

20.33 

Intermittency 

/(%) 

70.76 

70.70 

70.17 

66.49 

59.54 

42.30 

Relative  intensity 

7r 

4.90 

4.80 

4.50 

3.82 

2.95 

2.14 

